module equations
use technical
use grid
use mod_2dflu
use fluid
use mod_part_interp
use force
implicit none
public
save
!--------------------------
namelist /equations_param/ &
      lfluid,lparticle,lmagnetic,lpscalar,lpoly
contains
!***********************************************************!
subroutine initialize_equations
!
  tvar=0;naux=0
  call announce('initialize equations:','The details of equations being solved are as follows:')
  if(lfluid) then
! only variable that changes with time is vorticity
    tvar=tvar+1
    iomega=tvar
! auxiliary variables used are (ux,uy and psi) 
    naux=tvar+3
    call log_message('Solving for fluid:')
    call log_message('time-dependent variable: vorticity, tvar contribution = 1')
    call log_message('Auxilliary variable: ux,uy,psi, naux contribution = 4')
  else
  endif
!
! If there are more equations, similar analysis is done above to get the total
! number for tvar and naux.  At this stage we should know their values. 
!
  write(*,*)'log:','tvar=',tvar,'naux=',naux
  allocate(f(n1+2,n2,tvar),df(n1+2,n2,naux))
!
! Now call individual modules which sets up pointers to this array
!
  call point_fluid()
!
!
!
!
endsubroutine initialize_equations
!***********************************************************!
subroutine get_df(istep)
!
! Constructs the RHS of the equations we want to solve using the spcetral
! method. 
!
  integer :: istep
  integer :: i1,i2,ireal,iimag,k1,k2
!
! First ensure that the time-evolving part of df is set to zero.
!
  df(:,:,1:tvar) = 0.
!
!Calculation of df is done in two parts. In the first part all the quantities are
! in Fourier space. And corresponding Fourier space diagnostics are also 
! calculated. At the end of this part all quantities are transformed to 
! real space. 
!
! The following loop calculates all quantities local in Fourier space
!
  do i2 = 1,n2
     k2 = (i2-1) - n2*(i2/(n1hf+1))	
     do i1 =1,n1hf
        k1 = i1 -1
        ireal = 2*i1-1
        iimag=2*i1
        if(lfluid) & 
             call get_df_fluid_part1(istep,ireal,iimag,i2,k1,k2)
        if(lforcing_determinisic_fourier) & 
             call forcing_deterministic_fourier(ireal,iimag,i2,k1,k2)
     enddo
  enddo
! 
!Now FFT to obtain real space quantities
!
  if(lfluid) call fft_fluid_fourier_to_real ()
!
!
! The velocities are in real space. 
!
  if (lparticle) call dxp_dt(istep)
  if (lpoly) call dpoly_dt(istep)
!
! Calculations in real space
! (Once we include magnetic and passive scalar part, hooks that uses real
! space velocity should go here.)
  do i1=1,n1;do i2=1,n2
     if (lfluid) call df_fluid_real_space (istep,i1,i2)
  enddo;enddo
!
! Now FFT from Fourier to Real space
!
  if (lfluid) call fft_fluid_real_to_fourier ()
!
! Now do the second part of the get_df subroutine
!
  do i2 = 1,n2
     k2 = (i2-1) - n2*(i2/(n1hf+1))	
     do i1 =1,n1hf
        k1 = i1 -1
        ireal = 2*i1-1
        iimag=2*i1
        if(lfluid) &
             call get_df_fluid_part2(istep,ireal,iimag,i2,k1,k2)
     enddo
  enddo
!
!  
endsubroutine get_df
!***********************************************************!
subroutine read_equationsparam(fpointer)
  integer :: fpointer,ierr
  read(fpointer,NML=equations_param, IOSTAT=ierr)
  if(ierr/=0) call fatal_error('equations','problem with reading namelist')
endsubroutine read_equationsparam
!***********************************************************!
subroutine cleanup_equations
!
!cleans up allocated variables etc. 
!
  deallocate(f,df)
endsubroutine cleanup_equations
!***********************************************************!
endmodule equations
